The statistical mechanics of traveling salesman type problems 
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We study the finite temperature statistical mechanics of Hamiltonian paths between a set of N 
quenched randomly distributed points in a finite domain T>. The energy of the path is a function of 
the distance between neighboring points on the path, an example is the traveling salesman problem 
where the energy is the total distance between neighboring points on the path. We show how the 
system can be analyzed in the limit of large N without using the replica method. 
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In this Letter we consider a system where N points 
with quenched positions ar^ , x^ ■ ■ ■ Xjf are indepen- 
dently distributed on a finite domain T> with a proba- 
bility density function p q (x). In general, the domain T> 

is multidimensional and the points xf are vectors in 
the corresponding Euclidean space. Inside the domain T> 
we consider a polymer chain composed of N monomers 
whose positions are denoted by x\, x%, ■ ■ ■ , xn- Each 
monomer a;, is attached to one of the quenched sites 
x[ q ^ and only one monomer can be attached to each 
site. The state of the polymer is described by a permu- 
tation a G S^v where Sat is the group of permutations of 
N objects. The position of the i th monomer may thus be 
written as Xi — 3^2) ■ The Hamiltonian for the system is 



given by 
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Here V is the interaction between neighboring monomers 
on the polymer chain. For convenience the chain is taken 
to be closed, thus we take the periodic boundary condi- 
tion xq — xn- 

A physical realization of this system is one where the 
x^ are impurities where the monomers of a polymer loop 
are pinned. The potential V represents effective interac- 
tion between neighboring monomers on the chain. For 
instance V(x) = Ax 2 /2 corresponds to the Rouse model 
of a polymer chain[jj. Another application of this model 
occurs in combinatorial optimization: if one takes V(x) 
to be the norm, or distance, of the vector x then H(cr) is 
the total distance covered by a path which visits each site 
x[ q ^ exactly once. The problem of finding a* which min- 
imizes H(a) is known as the traveling salesman problem 
(TSP) |2(. If one takes V(x) = \x\ (the Euclidean norm) 
then the combinatorial optimization problem is called the 
Euclidean TSP [3(. The choice V(x) — —\x\ means that 
one is looking for the longest path, this is the so called 
Maximal TSP or taxicab rip-off 0. The analogy with 
the physical polymer system is particularly useful. First 
one may use the techniques of simulated annealing Q| 



to search for the optimal path a* by introducing a tem- 
perature in the system and then incorporating a dynam- 
ics respecting detailed balance. The optimal path corre- 
sponds to the ground state or zero-temperature energy 
Eqs of H . By slowly cooling the system one may find 
this state, though in the presence of many local minima 
or metastable states this procedure may not be efficient. 
Secondly, and this is the approach taken here, one may 
try to determine the full temperature dependence of the 
free energy averaged over the disorder ensemble and then 
take the zero temperature limit to determine the average 
value Eqs [EIEQ Depending on the form of V, Eqs 
may not be extensive in N. In such cases, for example for 
the Euclidean TSP, the temperature needs to be rescaled 
with TV in order to have an extensive ground state energy 
IH . In this letter we study the problem with a fixed 
domain size and no rescaling of the temperature with N. 
The TSP is often studied in this form and in the thermo- 
dynamic limit the ground state energy per site is strictly 
zero. 

The main technical problem associated with the sec- 
ond approach above is to perform the average over the 
quenched disorder. The replica method and cavity meth- 
ods have been used previously to perform the quenched 
average for the TSP g H H 111 Gj , notably in the ran- 
dom link version of the problem where the distances be- 
tween sites i and j are assumed to be independent. This 
'random link' assumption considerably facilitates carry- 
ing out the disorder average within the replica/cavity 
formalism. However, in the Euclidean TSP problem the 
distances between the sites are evidently correlated (for 
instance the triangle inequality must be respected). In 
this Letter we present an exact approach that (i) does 
not require the use of the replica method, (ii) fully takes 
into account the correlations between the distances and 
(iii) moreover provides us with exact asymptotic results 
at all temperatures. 

The canonical partition function is given by the sum 
over all permutations 



Zn = 



^ exp(-/?if(a)). 
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Since the number of permutations grows as N\, the en- 
tropy is of order iVlniV and we insert a factor of l/N\ to 
absorb it. To simplify our formulas we take unit domain 
size. 

We define the density of quenched sites on V as 



N 



Pq{ X ) = Jj^2 S ( X - X i q) ) 



(3) 
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The partition function Zjq only depends on p q {x). If X{ 
is the site visited by the monomer i then the partition 
function of the permutation problem can be written up 
to constant factors as 



Z N — 



i r N 



exp 




The delta function constraint above ensures that the 
monomers Xi have the same density as the quenched 
points and thus visit only the quenched points and with 
the correct degeneracy. Using a Fourier representation of 
the functional constraint gives, up to temperature inde- 
pendent constants, 



d\pt\ exp N / dxfi(x)p q (x) I Z^. (5) 



The object Zn is the annealed partition function for a 
free ring polymer whose monomers can attain any point 
in T> but with an x dependent chemical potential. It is 
defined as: 
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The partition function Zn may be evaluated by opera- 
tor techniques: Zj^ = Tr T N where T is the symmetric 
operator 



T(x, y) = exp 



(j,(x) p{y) 



- PV{x - y) 



The full partition function Eq. JSJ can be evaluated by 
the saddle point method in the limit where N — > oo keep- 
ing the size of T> fixed. The saddle point equation is 
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-(Y J 5{x~x l ))= Pa (x). (8) 



where the above expectation, is for the system with 
partition function Z^ defined in Eq. JHJ), so p a (x) is 
the, annealed, density of points (monomers) for the free 
ring polymer. Physically this approach can be thought 
of as choosing a site dependent chemical potential \x 
which fixes the density of the annealed calculation to 



be the same as that of the quenched one, i.e. so that 
Pq(x) = p a (x). The approach has some similarity with 
the constrained annealing approximation [l4L lig, but 
fixes the whole distribution rather than individual mo- 
ments. 

The ground state eigenfunction, corresponding to the 
maximal eigenvalue Aq of T obeys 



/o (9) 0*) = V / dyT(x,y)C(y) 



(9) 



and for large N the annealed density of points in T> is 
given by 



Pa(x) 



ain(A ) 



(10) 



Substituting this into the saddle point equation and writ- 
ing exp(-^i) = y/p q (x)/s\ (x) wc find that s Xo (x) 
obeys 



s\ Q (x) = A 1 dy eiq>(-f3V(x - y)) 



Pg(y) 
s\ (y)' 



(ii) 



Considering the case of a uniform distribution of the 
points, x\ q \ and substituting back into the action we 
obtain 



N 



N 



dx ln(sA (x)) +ln(Ao)+terms indep. of /3. 

(12) 

From Eq. (|11|) we see that there is a family of solutions 
{s\ (x), Ao} which are related by s\ = a}/ 2 s a \ , for a > 
and in addition these solutions all have the same action. 
This apparent zero mode is an artifact introduced by the 
fact that the constraint N = N J dx p q = J dx ^ 4 6(x — 
Xi) is automatically satisfied. Thus we may chose Ao = 1, 
leading to our final result for the average energy per site 



- = '"dp 



dx In (s(x)) 



dx dy fa^M , (13) 
s(x)s[y) 



(7) where s obeys 



s(x) 



dy 



exp(-(3V(x -y)) 



(14) 



In general Eq. i|14|) can be solved by an iterative nu- 
merical procedure. The average energies obtained in this 
way for the one dimensional TSP on the unit interval 
[0, 1] and the conventional two dimensional TSP on the 
square domain [0,1] 2 are shown as continuous lines in 
Fig. (1). To test these predictions we have carried out 
Monte Carlo simulations of this TSP for system sizes of 
N = 2000 and compared the average energy measured af- 
ter equilibrating the system over 10 6 Monte Carlo steps 
and measuring the average energy over a subsequent 10 6 
Monte Carlo steps. The basic move in the dynamics was a 
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FIG. 1: Theoretical prediction for the average energy e for 
one and two dimensional TSP as a function of (solid line) 
compared with the Monte Carlo simulations (solid circles). 
Negative (3 corresponds to the maximal TSP problem. Er- 
ror bars based on 20 realizations of the quenched points are 
smaller than the symbol sizes. 



transposition of a pair of points in the permutation, and 
the acceptance criterion was the Metropolis rule. The 
simulation points are also shown as solid circles on Fig. 
(1) and we see that for all temperatures the agreement 
with the theoretical prediction is excellent. 

For the one dimensional ring and two dimensional torus 
or sphere, where the Euclidean distance is taken to be 
the shortest way round the domain, Eq. I|14f> admits a 
constant solution. This indicates that in these cases, the 
annealed approximation is exact without the need for any 
chemical potential to force the density of the annealed 
points to be uniform. In the case of the ring we find 



1 1 

6 ~ P ~ 2(eP/ 2 - 1)' 



(15) 



for both positive and negative /3. This result, and the 
corresponding one for the torus and sphere are confirmed 
by Monte Carlo simulations, and moreover support the 
correctness of the iterative procedure used to numerically 
solve Eq. Ifn jl [Ta j . 

In one dimension the Euclidean TSP amounts to a 
sorting problem, but even when the quenched points are 
regularly spaced, no simple expression for the partition 
function is known [Tfij . For the unit interval V = [0, 1] 
we have obtained the low temperature expansion 



V(x) = — ln(|x|), the latter was only considered at posi- 
tive temperature as it is ill defined at negative tempera- 
ture. In all these cases, the comparison with simulation 
results confirms the predictions of our method. 

From an optimization point of view, besides the al- 
gorithm for discovering it, the most interesting quantity 
is the length or energy of the optimum path i.e. Eqs] 
thus we now focus on the zero temperature limit. First, 
let us recall that our treatment computes an extensive 
ground state energy. On the other hand, we expect that 
for the Euclidean TSP in d dimensions, the length of the 
shortest path should scale as N 1 ~ 1 / d 3]. In this Letter, 
rather than attempt to pursue finite size corrections to 
our formalism we study models which do have an 
extensive ground state energy. This is the case in the 
maximal TSP problem and for problems where the po- 
tential V(x) between neighboring monomers is repulsive. 
We reformulate Eq. (|13|) to work directly at zero tem- 
perature by writing s(x) = exp (— f3w(x)). Then in the 
limit (3 —> oo we find 



w(x) = min^o,!) {V(x - y) - w(y)} 



(17) 



The ground state energy is then given according to Eq. 
(E2J as 



(18) 



ecs = 2 / dx w(x). 



The Eq. (|17J) seems in general quite difficult to solve, 
however we can make progress in some specific cases in 
one dimension. First there is the obvious case of what 
happens when the minimum value of V(x) occurs at x = 
0, as is the case for the positive temperature TSP. Clearly 
£ gs = V(0) here as the optimal solution connects all the 
sites in ascending (or descending) order and the distance 
between each site ~ 1/N — > 0. In this case we see that 
Eq. (|17|l clearly has a solution w(x) = V(0)/2 and thus 
Eq. {18} implies indeed that cgs = V(0). 

When V is a monotonically decreasing function (so the 
monomers prefer to be as far from each other as possible) 
a possible strategy for obtaining the ground state energy 
is a greedy algorithm where one starts at the leftmost 
point and goes to the rightmost point, then to the second 
leftmost point and so on. This leads to an energy of 



dx V(x) 



(19) 



Obviously the solution w(x) to Eq. 1171) should be sym- 
metric about x = i thus w(x) = u(\x — ^|) and u obeys 
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(16) 



which agrees with our numerical results. We have also 
considered the one-dimensional potentials V(x) — x 2 and 



u(x) = mm ye( _ i _ i) {V(x - y) - u{y)} . (20) 

Inspired by the greedy algorithm we postulate that the 
minimization on the right-hand side of Eq. I|20l) is 
achieved by y — —x. This leads to a solution of Eq. 



P7| u(a 



hV(2x), which is valid when V"{2x) < 0, 



■i.e. 



when V is concave, for all x £ [0, 1]. The energy of this 



solution is given by Eq. l|19fl and is thus attained by the 
greedy algorithm. We cannot guarantee that this solu- 
tion is unique but it is in agreement with the simulations 
for the potentials V(x) — — \x\ and — x 2 |l3j . 

Another ansatz for solving Eq. (|20() is u(x) = a\x\ + b. 
The Eq. |JU is solved by this form when V"(l/2) > 
(and thus does not work for concave potentials) and 
one finds a = V'{\) and b = \ (V(§) - \V{\)). The 
ground state energy predicted by this solution is cgs = 
V{\) A potential where the above solution is possible 
is V(x) = — ln(|a;|). Numerical solution of Eq. (|14|> at 
low temperatures converges to the solution found above. 
The predicted value of the ground state energy is ecs — 
ln(2), this value is compatible with the simulations for 
this potential This energy is achieved by making 

jumps from site to site where the jump size is a random 
variable A very close to 1/2 i.e. if the current position 
is x < 1/2 one jumps to x + A, otherwise one jumps to 
x — A. This halfjump-algorithm will generate an energy 
per site of cha — V"(l/2) and it will also generate the 
required uniform quenched distribution of points on [0, 1]. 
To summarize, the greedy and half-jump algorithms will 
achieve energies per site: 

e GA = j dx V(x); e HA = V(±). (21) 

Clearly when V"(\x\) > everywhere in [0,1], Jensen's 
inequality tells us that (V(X)) > V((X)) for X dis- 
tributed on [0, 1]; when this distribution is uniform this 
implies that ea A > (ha and hence the half-jump algo- 
rithm is the most efficient. In the case where the potential 
is concave the greedy algorithm is the most efficient. We 
note that the case of the maximal TSP is an intermedi- 
ate case where V"(x) = and in this case to A = (ha 
and the forms of u(x) in these two cases coincide. If 
V(x) (for x > 0) is such that it has a single minimum 
at x = x* and where x* < 1/2, then Eq. l(T7|l has a 
solution w(x) = V(x*)/2 for all x and the ground state 
energy is eqs = V{x*). An algorithm which achieves this 
is the x*-jump algorithm, note that this works because 
when x* < 1/2 the algorithm be applied from any start- 



ing position. When x* > 1/2, w(x) — V(x*)/2 is clearly 
not a solution for all x as it fails in the neighborhood of 
x = 1/2. We note that when x* < 1/2 the jumping algo- 
rithm is in fact a greedy algorithm as it always chooses 
the near optimal jump size and this jump size is indepen- 
dent of its current position. When the optimal jump size 
x* > 1/2 the strategy clearly does not work for the rea- 
sons discussed above. In this case the jump sizes are no 
longer concentrated around some typical value and the 
next jump size depends on the current position. 

We have discussed the statistical mechanics of models 
whose phase space is the set of permutations of N objects 
characterized by quenched positions/sites. The Hamilto- 
nians are functions of the neighboring elements in the 
sequence, and thus a given sequence can be interpreted 
as the energy of a polymer or random walk which visits 
each site once. We showed how the quenched calculation 
could be carried out and confirmed its predictions with 
Monte Carlo simulations. Physically the method corre- 
sponds to imposing a fictitious site dependent chemical 
potential on the distribution of the set of dynamical vari- 
ables Xi in the presence of the original interaction Hamil- 
tonian. This chemical potential is then chosen to ensure 
that the annealed density of these dynamical x\ is the 
same as the desired distribution of the quenched random 
variables x^ . The method introduced works in the ther- 
modynamic limit (corresponding to high density where 
the size of the domain is held constant) for any quenched 
distribution p q (x), in any dimension and for any inter- 
action potential V(x). For the problem of a directed 
polymer in dimensions greater than two a finite temper- 
ature phase transition is known to occur it would 
be interesting to see if this phase transition shows up in 
the Euclidean TSP in d > 2. Finally, the idea of treat- 
ing quenched variables as effectively annealed variables 
and then adjusting their Boltzmann weight in order to 
recover, self consistently, the original quenched distribu- 
tion may prove useful either as an exact or approximate 
method, as is the case of Morita's approach 0,13) m 
other problems involving quenched disorder. 
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